/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Jurgen Van Gael
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;
using dnAnalytics.Math;
using dnAnalytics.Properties;

namespace dnAnalytics.Statistics.Distributions
{
    /// <summary>
    /// The negative binomial is a distribution over the natural numbers with two parameters r,p. For the special
    /// case that r is an integer one can interpret the distribution as the number of tails before the r'th head
    /// when the probability of head is p.
    /// </summary>
    /// <remarks>The distribution will use the <see cref="System.Random"/> by default. 
    /// Users can set the random number generator by using the <see cref="RandomNumberGenerator"/> property.
    /// 
    /// The statistics classes will check all the incoming parameters whether they are in the allowed
    /// range. This might involve heavy computation. Optionally, by setting Control.CheckDistributionParameters
    /// to false, all parameter checks can be turned off.</remarks>
    public class NegativeBinomial : IDiscreteDistribution
    {
        private readonly double mP;
        private readonly double mR;

        /// <summary>
        /// Construct a new negative binomial distribution.
        /// </summary>
        /// <param name="r">The r parameter.</param>
        /// <param name="p">The p parameter.</param>
        public NegativeBinomial(double r, double p)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(r, p);
            }

            mR = r;
            mP = p;
        }

        /// <summary>
        /// The r parameter.
        /// </summary>
        public double R
        {
            get { return mR; }
        }

        /// <summary>
        /// The p parameter.
        /// </summary>
        public double P
        {
            get { return mP; }
        }

        /// <summary>
        /// A string representation of the distribution.
        /// </summary>
        public override string ToString()
        {
            return "NegativeBinomial(R = " + mR + ", P = " + mP + ")";
        }

        #region IDistribution Members

        /// <summary>
        /// The mean of the distribution.
        /// </summary>
        public double Mean
        {
            get { return mR*(1.0 - mP)/mP; }
        }

        /// <summary>
        /// The standard deviation of the distribution.
        /// </summary>
        public double StdDev
        {
            get { return System.Math.Sqrt(mR*(1.0 - mP))/mP; }
        }

        /// <summary>
        /// The variance of the distribution.
        /// </summary>
        public double Variance
        {
            get { return mR*(1.0 - mP)/(mP*mP); }
        }

        /// <summary>
        /// The entropy of the distribution.
        /// </summary>
        /// <value></value>
        public double Entropy
        {
            get { throw new Exception("Not implemented yet."); }
        }

        /// <summary>
        /// Gets or sets the random number generator.
        /// </summary>
        /// <value>The random number generator used to generate a random sample.</value>
        public System.Random RandomNumberGenerator { get; set; }

        #endregion

        #region IDiscreteDistribution Members

        /// <summary>
        /// The mode of the distribution.
        /// </summary>
        public int Mode
        {
            get { return (mR > 1.0 ? (int) System.Math.Floor((mR - 1.0)*(1.0 - mP)/mP) : 0); }
        }

        /// <summary>
        /// The median of the distribution.
        /// </summary>
        /// <value></value>
        public int Median
        {
            get { throw new Exception("Not implemented yet."); }
        }

        /// <summary>
        /// Computes the probability of a specific value.
        /// </summary>
        public double Probability(int val)
        {
            double ln = SpecialFunctions.GammaLn(mR + val)
                        - SpecialFunctions.GammaLn(mR)
                        - SpecialFunctions.GammaLn(val + 1.0)
                        + mR*System.Math.Log(mP)
                        + val*System.Math.Log(1.0 - mP);
            return System.Math.Exp(ln);
        }

        /// <summary>
        /// Samples a negative binomial distributed random variable.
        /// </summary>
        public int Sample()
        {
            return DoSample(RandomNumberGenerator, mR, mP);
        }

        /// <summary>
        /// Samples a negative binomial distributed random variable.
        /// </summary>
        /// <param name="n">The number of variables needed.</param>
        public int[] Sample(int n)
        {
            return DoSample(RandomNumberGenerator, n, mR, mP);
        }

        #endregion

        /// <summary>
        /// The skewness of the distribution.
        /// </summary>
        public double Skewnewss
        {
            get { return (2.0 - mP)/System.Math.Sqrt(mR*(1 - mP)); }
        }

        /// <summary>
        /// Check the parameters of the distribution.
        /// </summary>
        /// <exception cref="ArgumentOutOfRangeException">If r is nonpositive or p is not in [0,1].</exception>
        private static void CheckParameters(double r, double p)
        {
            if (r < 0.0)
            {
                throw new ArgumentOutOfRangeException("r", Resources.ZeroInfinityRange);
            }

            if (p < 0.0 || p > 1.0)
            {
                throw new ArgumentOutOfRangeException("p", Resources.ZeroOneRange);
            }
        }

        /// <summary>
        /// Samples a negative binomial distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="r">The mean of the Laplace distribution.</param>
        /// <param name="p">The scale of the Laplace distribution.</param>
        public static int Sample(System.Random rnd, double r, double p)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(r, p);
            }

            return DoSample(rnd, r, p);
        }

        private static int DoSample(System.Random rnd, double r, double p)
        {
            double lambda = Gamma.Sample(rnd, r, p);
            return Poisson.Sample(rnd, lambda);
        }

        /// <summary>
        /// Samples a negative binomial distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of variables needed.</param>
        /// <param name="r">The r parameter.</param>
        /// <param name="p">The p parameter.</param>
        public static int[] Sample(System.Random rnd, int n, double r, double p)
        {
            if (Control.CheckDistributionParameters)
            {
                CheckParameters(r, p);
            }

            return DoSample(rnd, n, r, p);
        }

        /// <summary>
        /// Samples a negative binomial distributed random variable.
        /// </summary>
        /// <param name="rnd">The random number generator to use.</param>
        /// <param name="n">The number of variables needed.</param>
        /// <param name="r">The r parameter.</param>
        /// <param name="p">The p parameter.</param>
        private static int[] DoSample(System.Random rnd, int n, double r, double p)
        {
            int[] arr = new int[n];
            for (int i = 0; i < n; i++)
            {
                arr[i] = DoSample(rnd, r, p);
            }
            return arr;
        }
    }
}